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Abstract 


Research in psychology and neuroscience has successfully modeled decision making as a pro¬ 
cess of noisy evidence accumulation to a decision bound. While there are several variants 
and implementations of this idea, the majority of these models make use of a noisy accumu¬ 
lation between two absorbing boundaries. A common assumption of these models is that 
decision parameters, e.g., the rate of accumulation (drift rate), remain fixed over the course 
of a decision, allowing the derivation of analytic formulas for the probabilities of hitting the 
upper or lower decision threshold, and the mean decision time. There is reason to believe, 
however, that many types of behavior would be better described by a model in which the 
parameters were allowed to vary over the course of the decision process. 

In this paper, we use martingale theory to derive formulas for the mean decision time, 
hitting probabilities, and first passage time (FPT) densities of a Wiener process with time- 
varying drift between two time-varying absorbing boundaries. This model was first studied 


by Ratcliff (1980) in the two-stage form, and here we consider the same model for an arbi¬ 


trary number of stages (i.e. intervals of time during which parameters are constant). Our 
calculations enable direct computation of mean decision times and hitting probabilities for 
the associated multistage process. We also provide a review of how martingale theory may 
be used to analyze similar models employing Wiener processes by re-deriving some classical 
results. In concert with a variety of numerical tools already available, the current derivations 
should encourage mathematical analysis of more complex models of decision making with 
time-varying evidence. 
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1. Introduction 


Continuous time stochastic processes modeling a particle’s diffusion (with drift) towards 
one of two absorbing boundaries have been used in a wide variety of applications including 
statistical physics (Farkas and Fulop, 2001), finance (|Lin 1998), economics (Webb, 2015), 


and health science (Horrocks and Thompson 2004). Varieties of such models have also been 


applied extensively within psychology and neuroscience to describe both the behavior and 
neural activity associated with decision processes involved in perception, memory, attention, 


and cognitive control ( 
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, 1992 I 

datediff and Rouder, 1998 d 

tatcliff and McKoon 

2008] 

Simen et al 

2009] Gold anc 

Shadlen, 2001, 2007] E 

Irunton et a 
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2009] 
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2010). 


In these stochastic accumulation decision models, the state variable x(t) is thought to 
represent the amount of accumulated noisy evidence at time t for decisions represented by 
the two absorbing boundaries, that we refer to as the upper (+) and lower (—) thresholds 
(boundaries). The evidence x(t) evolves in time according to a biased random walk with 
Gaussian increments, which may be written as dx(t) ~ Normal(/i dt, a 2 dt ), and a decision 
is said to be made at the random time r, the smallest time t for which x(t) hits either 
the upper threshold (x(t) = +C) or the lower threshold (x(r) = — (), also known as the 
first passage time (FPT). The resulting decision dynamics are thus described by the FPT 
of the underlying model. In studying these processes one is often interested in relating the 
mean decision time and the probability of hitting a certain threshold (e.g. the probability of 
making a certain decision) to empirical data. For example, these metrics can offer valuable 
insight into how actions and cognitive processes might maximize reward rate, which is a 


simple function of the FPT properties (Bogacz et ah, 2006). 


However, not all decisions can be properly modeled if parameters are fixed throughout the 
duration of the decision process. Certain contexts can be better described by a model whose 
parameters change with time. In this article we analyze the time-dependent version of the 
Wiener process with drift between two absorbing boundaries, building on recent work that 


is focused on similar time-varying random walk models (Hubner et ah, 2010 Diederich and 


Oswald, 2014). After reviewing how martingale theory can be used to analyze and re-derive 


the classical FPT results for the time independent case, we calculate results for the time- 
dependent case. The main theoretical results are presented in §5.2[ where we provide closed 
form expressions for threshold-hitting probabilities and expected decision times. In | Appendix 


[P| we also describe how our methods can be applied to the more general Ornstein-Uhlcnbeck 
(O-U) processes, which are similar to the Wiener diffusion processes albeit with an additional 
“leak” term. We conclude with a summary of the results and a discussion of how the present 
work interfaces with other similar analyses of time-varying random walk models. 


2. Notation and terminology 

Here we introduce the notation and terminology for describing the model we analyze, 
which is a Wiener process with (time-dependent) piecewise constant parameters. This sim¬ 


ple stochastic model, and others close to it, have been studied before (Ratcliff, 1980 Heath 
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1992 

Smith 

2000 

Diederich and Busemeyer 

2003, 2006; Diederich and Oswald 

2014; 

Bo- 

gacz et al. 

2006 

Wagenmakers et al. 

2007]), although the reader should note that our 


parameterization differs from that of some previous studies. Before describing the model, 
we first review the case where parameters are unchanging with time. In order to easily 
discuss this simpler model alongside the main time-dependent model analyzed in (|5j we 
use the terms single-stage model and multistage model , respectively. Readers familiar with 
the popular Diffusion Decision Model of Ratcliff and McKoon (20083 should be aware that 
parameters in our model do not vary randomly trial-to-trial. Readers familiar with 
(2006) should be aware that the single-stage model Q is equivalent to what 


Bogacz 


et al. 


et al. 


Bogacz 


(2006) call the “pure” Drift Diffusion Model. 


2.1. The single-stage model/process with constant parameters 
Consider the stochastic differential equation (SDE): 

dx(t) =/a dt + crdW(t), x(to) = xq, 


( 1 ) 


where parameters /i and a are constants referred to as the drift and diffusion rates, re¬ 
spectively; Xq is the initial condition (the initial evidence or starting point) of the decision 
process; and adW{t) are independent Wiener increments with variance a 2 dt. This simple 
stochastic model has successfully modeled the evolution of evidence between two decisions 


during a two-alternative forced choice task (Ratcliff and Rouder, 1998 Bogacz et al. 2006 


Ratcliff et al. 2016), so that (|TJ) can be interpreted as modeling a decision process in which 


an agent is integrating noisy evidence until sufficient evidence is gathered in favor of one of 
the two alternatives. 

A decision is made when the evidence x(t ) crosses one of the two symmetric decision 
thresholds ±£ for the first time, also referred to as its first passage time (FPT). In other 
words, a decision occurs the instance x(t) crosses and is absorbed by one of the two bound¬ 
aries. The two boundaries each correspond to one of the two possible decisions for the task. 
We will refer to the absorbing thresholds at +£ and — ( as the upper and lower decision 
boundaries/thresholds. To contrast with the next section, we will sometimes refer to this 
model with time-invariant parameters as the single-stage model or single-stage process. We 


again note that the parameterization used here differs from that employed by others (Smith 


2000 Ratcliff and Smith, 2004 Navarro and Fuss, 2009), although the underlying model is 


equivalent. Our formulation, compared to some others, does not include a parameter for 
“non-decision time” or “timeout” for a given trial. Such a term could be incorporated by 
shifting the entire reaction time distribution - it has no effect on any of our analyses. 

2.2. Time-dependent, piecewise constant parameters 

The assumption that model parameters remain constant throughout the decision process 
is unlikely to hold in many situations. For example, the quality of evidence may not be 
stationary (i.e., the drift rate and diffusion rate are not a constant with respect to time) 
or decision urgency may require thresholds to collapse in order to force a decision by some 


2012 ). 


deadline (Cisek et al. 2009 Mormann et al. 2010; Zhang et al. 2014 Drugowitsch et al. 
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In order to analyze such situations, we focus our present study on a two-stage model 
originally analyzed by Ratcliff (1980), which we generalize to an arbitrary number of stages. 
We refer to this slightly generalized model as a multistage model or multistage process , to 
distinguish it from ([Tj) above. The multistage model allows for the drift rate, diffusion rate, 
and thresholds to be piecewise constant functions of time. 

To fully describe the multistage model, we first partition the set of non-negative real 
numbers (i.e. time axis) into n intervals (or stages ) [L_i,tj], i e {1,... ,n} with t 0 = 0 and 
t n = Too. We then assume that the drift rate, the diffusion rate, and the decision thresholds 
are constant within each interval, but that their values change between intervals. Evidence 
integration is thus modeled by 


dx(t) = n(t) dt T a(t)dW(t), x(to) = xq, 


( 2 ) 


where 


A i(t) = /h, for U -1 <t <U , 
cr{t) = cTj, for T_! < t < ti , 


for each i e {1, 


, n). The above assumptions are identical to the assumptions in Diederich 
and Qswald| ([2014 2016). If n — 1, the multistage model reduces to the single-stage 
model ([!]). For expository clarity, we begin by assuming the decision thresholds are fixed at 
±£, and in qOwe generalize to time-varying (piecewise constant) thresholds. Let r be the 


( 3 ) 


first passage (decision) time for the multistage model. 

We will frequently refer to the i-th stage of (|2]), which for t > ti- 1 , is written as 

dx(t) = a kdt T < TidW(t ), x(U-i) = W_i, decision thresholds ± (, 


where the initial condition X Wi is a random variable dehned as x{ti-\) conditioned on there 
being no decision (threshold-crossing) before time More precisely, the density of Xi_i is 
the conditional distribution of x(T_i) given that r > ti- 1 - Thus, the random variable W-i 
corresponds to realizations of the multistage model that remain within the thresholds ±£ 
until time U-i. For the first stage, X 0 could be a point mass centered at ^o, or it may be 
a random variable capturing the variability in starting points (Ratcliff and Rouder, 1998). 
The key difference between (|3]) and the single-stage process (JT|) is that the initial condition 
Xi-i is a random variable whose distribution is determined by previous stages. 


3. Martingale theory applied to the single-stage model 

In this section, we give an introduction to the basic properties of martingales and the 
optional stopping theorem, which are the key mathematical tools used in calculating our 
main results in §[5] For readers who are less familiar with martingale methods, we first 
derive the mean decision time, hitting probabilities, and FPT densities for the single- and 
two-stage models. These analyses provide an alternate approach to deriving these classical 


results as compared to other non-mart ingale based approaches (Ratcliff, 1980 Diederich and 


Oswald, 2014). We discuss the differences between these various approaches in 
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3.1. Continuous time martingales 

Consider a continuous time stochastic process 77 (f), t > 0. Let rj(t i : ^ 2 ) denote the 
portion of r](t ) between times G and ^ 2 - A stochastic process M(f) is said to be a martingale 
with respect to r](t) if the following three conditions hold: 

1. M(t) is a function of r/(0 : t ) and does not depend on future values of 

2. K\M(t)\ < 00 

3. M[t ) satisfies the stationarity condition in expected value 

E[M(f 2 )|^(0 : fi)] = M(ti), for all 0 < ti < t 2 < 00 . (4) 

The first condition means that given the realized values r)(0 : t ), we should be able to 
compute M(t ) deterministically, so that M(t) does not depend on the future. The second 
condition is a regularity condition that ensures that M(t ) is sufficiently light-tailed and 
holds under several decision-making scenarios. The third condition is the most crucial - it 
enforces stationarity in expected value. This last condition can be interpreted as a “fair 
play” condition ensuring that chances of gaining and losing starting with value M{t \) at 
time t\ are the same, as in the classic example of a sequence of flips of a fair coin. When 
introducing a martingale, one often does not explicitly specify the process git), and in this 
case 77 (f) is assumed equal to M(t). 

Martingale theory is very broad and there are many different choices for M(f) and 77 
which are interesting. All of the calculations and results of this paper are constructed from 
two fundamental stochastic processes, which we now introduce. Let W(t) be the standard 
Wiener process and X{t) = (it + aW{t) be a single-stage Wiener process with drift rate /i 
and diffusion rate a (without boundaries). We now consider some martingales associated 
with these two stochastic processes: 

1. W[t) is a martingale. It is easy to verify that for 0 < ti < f 2 , E[kF(f 2 )|IF(fi)] = W{ti). 
Similarly, X{t) — /it is a martingale. 

2 . W(t) 2 3 — t is a martingale. Note that for 0 < t\ < t 2 , conditioned on W(t\), W{t 2 ) has 
a Gaussian distribution with mean W{t\) and variance ( t 2 — t±). Therefore, 

E[W(t 2 ) 2 \W(t 1 )] = W(t 1 ) 2 + (t 2 -t 1 ) => E[W(t 2 ) 2 -t 2 \W(t 1 )} = W(t 1 ) 2 -t v 

3. For any A 6 R, exp(AX(f) — A /it — \ 2 a 2 t/2) is a martingale. Note that for 0 < ti < t 2 , 
conditioned on X(ti), X(t 2 ) has a Gaussian distribution with mean X(ti) + /i(f 2 — G) 
and variance a 2 {t 2 — t\). Thus, 

E[exp(AX(f 2 ))|W(fi)] = exp (A(X(G) + n(t 2 - G)) + A 2 a 2 (t 2 - G)/2) 
E[exp(AA"(f 2 ) — A 7 A 2 — X 2 a 2 t 2 /2)\X(ti)\ = exp (AA"(ti) — A/iti — A 2 a 2 fi/ 2 ) . 


1 More precisely, M(t) is progressively measurable with respect to the sigma algebra generated by 77(0 : t). 
chap. 2), (|Karatzas and Shreve 


See (Doob 


1953 


1998 


chap. 1), or (Revuz and Yor 


2 M ( t) may be heavy tailed due to non-existence of second and higher moments. 


1999 


chap. 1). 
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For A = —2 fi/a 2 , this martingale reduces to exp(—2 g,X{t)/o 2 ) which is referred to as 
the exponential martingale. 

3.2. Stopping times and the optional sampling theorem 

The first passage time r is a random variable defined by r — inf {t > 0| x(t) i (— C, C)}- 
We are interested in computing conditional expectations and probability densities of r, which 
correspond to expected decision times and the corresponding distributions of response times. 

The key tool we borrow from the theory of martingales is a classic result known as Doob’s 
optional sampling theorem (also known as the optional stopping theorem), which we motivate 
and introduce here. To understand the optional sampling theorem, one must first recall that 
the expected value of a martingale M(t) computed over all realizations starting from M(0) 
is equal to the initial expectation of M(0). That is, martingales by definition must satisfy 
the following: 

E [M(t)] = E[M(0)], for each t > 0. 

One then wonders: Does a similar property extend to the random time r? More specifically, 
if we consider different realizations of r and compute averages of M(r) at these realized 
values, does this average, as the number of realizations grow large, converge to M(0)? The 
answer is affirmative if r is well behaved and is independent of the process M(t). Indeed, in 
this case 

E[M(t)|M( 0)] = E[E[M(f)|M(0),r = t]\ = E[M(0)], (5) 

where the outer expectation is with respect to r and the inner expectation is with respect 
to M(t). 

But what if r is not independent of M(t)l In these cases the situation is more subtle. 
Suppose r is bounded from above by t. Then, we can write 

dM{u) = M{ 0) + [ 1 (t > u)dM{u ), 

Ju= 0 

where l(-) is the indicator function that takes value 0 if its argument is false and 1 otherwise, 
and dM{u) is the (random) increment in M{u ) at time u. If we assume that the value 
1 (t > u ) can be deterministically computed based on the knowledge of M(0 : u ), we can 
then write 

E[M(t)} = E[M(0)] + / E[1(t > u)dM{u)] = E[M(0)] + / E[E[l(r > s)dM{u)\M{0 : u)]], 

J u =0 J u =0 

where the second equality follows from the law of total expectation)^] The first equality 
requires swapping of integral and expectation operators, which is allowed because i is fi¬ 
nite. Furthermore, 1 (t > u) is a deterministic function of M(0 : u) and thus, E[l(r > 


M(t) = M{ 0) + 


' u=0 


3 For an integrable random variables Y and an arbitrary random variable Z 1 E[Y] = £ , [£’[Y|.Z]]. Loosely 
speaking, the law of total expectation states that the expectation of a random variable can be computed by 
first computing the expectation conditional on another random variable, and then computing the expected 
value of the resulting expectation. 
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u)dM(s)\M(Q : u)} = ~E[dM(u)\M(Q : w)]l(r > u) — 0, where the last equality follows 
by definition of martingale. Consequently, for a random variable r and martingale M(t), 
E \M(t)] = E[M(0)] if (i) the event r > u is determined by M(0 : u), and (ii) r is bounded 
from above with probability one. A random variable satisfying the first condition is called a 
stopping time, and the above discussion is the optional sampling theorem which we formally 
state: 

The optional sampling theorem. Suppose M(t ), t > 0 is a martingale with respect 
to rj(t) and r is a bounded (with probability one) stopping time with respect to git), then 
E[M(r)] = E[M(0)]. 


Heuristic-ally, the optional sampling theorem states that different realizations of a martin¬ 
gale M (■ t ) stopped at random times average out to constitute a fair game. The crucial aspect 
is that the stationarity of the expected value holds even for random (stopping) times, includ¬ 
ing our first passage time r. As we will see in 1R3 this stationarity enables us to calculate 
analytic expressions for first passage time properties by finding appropriate martingales. 

A helpful example is to consider the standard Wiener process with initial position at 
X(0) = xq, and absorbing thresholds at ±£, with x$ between ±£. In this case X(t) is 
itself a martingale, the first passage time r is a stopping time, and the optional sampling 
theorem says that E[X(0)] = E[A(r)]. The expectation on the left hand side is simply the 
average of the initial distribution of X (0) which is the number xq. The right hand side 
is more interesting: X{r) is the random value of X(t) at the random decision time r, the 
instant X(t) crosses +£ or Thus X(r) attains one of two possible values, +C or — (, and 
E[X(t)] heuristically resembles an average over all +£’s and — £’s corresponding to sample 
paths starting at Xo and diffusing until they hit either +£ or —£ at r. The optional sampling 
theorem says that this average of +C’s and —£’s, over all such sample paths, ends up being 
equal to the number xq. 


3.3. Applications to the single-stage model 

The optional sampling theorem is a powerful mathematical tool for decision-making mod¬ 
els that associate decisions and decision times with a diffusion processes crossing a threshold. 


Ratcliff’s Diffusion Decision Model ( 

Ratcliff, 

1978 

Ratcliff and McKoon, 

2008 

), the leaky 

competing accumulator model ( 

Usher and McClelland, 

2001 

), and the EZ diffusion model 


(Wagenmakers et al. 2007), are popular examples of such models. The optional sampling 


theorem reduces the problem of computing analytic expressions for the statistics of the first 
passage times to identifying appropriate martingales. In this section we illustrate the flavor 
of such calculations for the single-stage model from (JTJ) in §2.1 Recall that the decision 
time r is defined by r = inf (t > 0|a;(t) ^ (—C, C)}- Throughout this section, we introduce 
6 = (n,a, Q to slightly condense the notation when desired. 

We first compute P[.x(r) = —£], the probability of hitting the lower threshold. First, 
we let s = p/a 2 , the ratio of the drift parameter to the squared diffusion parameter (i.e. 
signal to noise). Recall from 13.1 that for p ^ 0, exp(— 2sX) is a martingale. Applying the 
optional sampling theorem, we get 


exp(—2sx 0 ) = E[exp(—2sa;(r))] = P[x(r) = C]E[exp(—2s£)] + P[x(r) = —C]E[exp(2s£)]- 


7 




























Substituting P[x(r) = £] = 1 — P[x(r) = —£] and solving for P[x(r) = —C], we obtain a 
closed form expression 

tot / \ = = exp(-2sx 0 ) - exp(-2sQ 

[ { ^ exp(2<) - exp(-2<) ' 

Similarly, for fi — 0, we note that X(t) is a martingale. Applying optional sampling theorem, 
we get 

x 0 = E[X(r)] = CP[x(r) = C] - CP[z(r) = -(], 

and following the same argument we obtain P [x(t) = —C] = (C — Xq)/2(. In summary, we 
get 

{ exp(— 2sx 0 ) - exp(—2Q , f ^ Q 

exp(2<) - exp(-2<) ’ ’ ( 6 ) 

where P ± (x o, 6) is the probability of hitting the upper and the lower threshold, respectively. 

To compute the expected decision time E[r] for /i ^ 0, recall from 13.1 that X(t) — fit is 
a martingale. For g ^ 0, applying the optional sampling theorem yields 

x 0 = E[(C - /xr)]P[x(r) = C] + E[(~C ~ /xr)]P[x(r) = -(]• 


Solving for E[r], we get E[r] = (1 — 2 P[x(t) = — C])C/h- When fi — 0, recall from 13.1 that 
X(t) — a\t is a martingale, and the same argument as above yields E[r] = (£ 2 — Xq )/a 2 . In 
summary, the mean decision time mDT(xo,6 ) ) is given by 


(1 - 2 P (x 0 ,6»))C - x 0 


mDT(xo,6 ) ) = E[r] = 


c 2 -^ 






if n 7^ 0, 
if fJL — 0. 


(7) 


We also wish to hnd r’s Laplace transform, or moment generating function^ E[exp(ar)]. 
We remember from 13.1 that exp(AX(t) — A/it — A 2 n 2 t/2) is a martingale, and choose A so 
that the coefficient of t becomes —a, i.e., A solves the equation |cr 2 A 2 + /iA — a = 0. The 
two solutions to this equation are 

—// — \jn + 2 aa 2 
Ai — -x-, 


and Ao = 


+ vV + 2a<j2 


(7 (7 

Applying the optional sampling theorem, we obtain a pair of equations, 

exp(AiXo) = exp(AiC)E[exp(—ar)l(x(r) = £)] + exp(—AiC)E[exp(—ar)l(x(r) = —£)] 
exp(A 2 x 0 ) = exp(A 2 C)E[exp(-ar)l(x(r) = £)] + exp(-A 2 C)E[exp(-ar)l(x(r) = — £)], 


4 The moment-generating function (technically, the two-sided Laplace transform) of a random variable X 
is <f>x(ot) = E[exp(aA')], a function of a £ A C R. It is often of interest because it specifies the probability 
distribution of X, and can be used to obtain the moments of X. 















which we can solve simultaneously for E[exp(— oit)1{x{t) = ±C)] to obtain 

/M±C-^hsmi y {fa ” ) 'fF T E 


E[exp(— ar)\{x{r) = ±C)] = exp 


cr- 


sm 


, f 2£-y/2a<x 2 +/i 2 - 


Thus, the moment generating function of the decision time is 


E[exp(— ar)\ = exp 


MC-*ohsmh( «+" ) 'fF^ ) , , M (C + xo)Nsmh( <Vgypg±g ) 




-) 


sinh( 


2£-y/2a:er 2 + /i 2 - 


+exp 


(7 




sinh( 


2aa~ 


As a byproduct, we also get the Laplace transform of conditional decision times: 

cp (MC^o}) sinh (A±Wpy!±Z) 


E[exp(-ar)|x(r) = C] = 


E[exp(-ar)|x(r) = -(] = 


exr 


P+( Xo ,6) 


exp(- 




) sinh( 


(C~^o)\/2acr 2 +/j 2 . 


P-(x o,0) 


. , /2C\/2 acr 2 +fj, 2 x 

smh(- SL j-) 


( 8 ) 


( 9 ) 


( 10 ) 


The derivatives of the Laplace transform yield moments of decision time (see Srivastava 


et al. (2016) for detailed derivation of conditional and unconditional moments of decision 


time using Laplace transforms). Here, we focus on expressions for conditional expected 
(mean) decision times that are the derivative of the Laplace transform with respect to —ex 
computed at a = 0. The expected decision time conditioned on hitting the upper and lower 
boundaries are denoted by mDT + and mDT“, and may be computed by differentiating (|9]) 
and flip]): 


mDT + (xo, 0) = 

mDH(x 0 ,^) = 


mDT + (x 0 , 9) \ f coth(2<) - coth(s(C + x 0 )), if fi 0, 


P+(x 0: 6) 

mDT (x 0 ,6*) 
P-{x 0 ,6) 


if fi — 0; 


4C 2 _ (C±goy \ 

3cr 2 3 <j 2 > ’ 

f coth(2<) - coth(s(C - x 0 )), if fi ± 0, 


((— x o) : 

3<j 2 


). 


if n — 0, 


in) 

( 12 ) 


where mDT = E[rl(x(r) = ±£)] and l(-) is the indicator function. We again note, just as 
with P[x(r) = —£], mDT ± also depend on the underlying parameters /x, xo,C> an d a - 

We now compute r’s first passage time density /(£; xq, /i, a, (), i.e., the probability density 
function of the decision time. This amounts to calculating the inverse Laplace transform of 
(|8|. In this case, the inverse Laplace transform needs to be expressed as an infinite series 
(see 


Lin (1998) for a detailed derivation): 


f(t;x Q ,6) = [t < t] = exp ( exp(-s(C + x 0 ))ti(t-, 


C~xq 2C 

a ’ cr 


exp(s(C - x 0 ))i9[t] ^ (13) 


cr 
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where 9 = 
defined by 


(a, a, () , and i9(t;u,v) is a function (Borodin and Salminen, 2002, pp. 451) 


+oo 

#(t;u,v) = 

k =—oo 


v — u + 2 kv 
f 3 / 2 


exp 


( (v — u + 2 kv) 2 
V 2t 


u < v. 


Similarly, the first passage time density conditioned on a particular decision is given by 


[T < t\x(T ) = C] 
jF[t < t\x(T) = -C] 


f + (t]x 0 ,6) 
P + {x o,9) 
f-(t-,x 0 ,6) 
P-(xo,0) 


exp(-g+ s(C-x 0 )) Q / C + ^o 2(\ 
P+(x 0 ,e) V ’ a ’ J 

ex P(“S ~ s (C + ^o)) Q f ( ~ xo 2(\ 

P~{x 0 ,9) V' a ’ J’ 


(14) 

(15) 


where / ± (t;xo,6 l ) = [r < t & x(r) = ±^], i.e., f ± (t]Xo,9)dt is the probability of the 
event r G [t,t + dt ) and x(r) = ±£. Note that / defined in (13) is the sum of f + and f~. 
Alternate derivations for the hitting probabilities, mean decision times, and FPT densities 


may be found in the decision making literature (Ratcliff and Smith, 2004; Bogacz et ah, 2006 


Navarro and Fuss, 2009). It is worth noting that the infinite series solution for the FPT 


density given in (13) is equivalent to the small-time representations for the FPT analyzed 
in (Navarro and Fuss (2009) and Blurton et al. (2012)). For completeness, we provide the 


alternative expression for density in |Appendix A 


4. Analysis of the two-stage model 

In this section, we use the tools developed in 


I in order to analyze the two-stage process 


originally presented and analyzed in Ratcliff (1980). While our calculations lead to equiv¬ 
alent formulas for the first passage time densities, a martingale argument provides us with 
additional closed form expressions for the probability of hitting a particular threshold and 


expected decision times. Computations of these FPT statistics using the results of Ratcliff 


(1980) requires numerical integration of the FPT density, which our formulas now avoid. 
We may explicitly write the two-stage model as 


dx{t) = dt + a(t)dW(t), x(t 0 ) = xq, 


where 


p(t) 



for 0 < t < 1 1 , 
for t > t\, 


and a(t) 


0 i, for 0 < t < ti, 
u 2 , for t > t\. 


(16) 


As before, the decision time r = inf {A > 0|x(t) ^ ( — CO} is th e fi rs t passage time with 
respect to boundaries at ±(. Let 9 % = (//,, cr*, ti, C) and s* = /q/of, for i G {1,2}. 
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We interpret the first stage as a single-stage model with a deadline at t\. For a single- 
stage model with thresholds ±£ and a deadline ti, the joint density g ddln (x,r]Xo,ti,9i) of 
the evidence x(ti) and the event t >t\ is given by (Douady, 1999 Durrett 2010): 


/ dln (x, r; xq, h, 9i) = -^-P[a;(ti) < x & r > ti] 


dx 
l\x E 


-c,0] 




exp 


-fi\t + 2gi(x - £ 0 ) 


2 of 


x 


E 


exp 


n=—oo 


— (x — xo + 4n£) s 
2ofi 


exp 


— (2£ — x — x 0 + 4nC) 5 
2aft 


(17) 


Here superscript “ddln” refers to the deadline. g ddlu may then be used to determine the 
FPT distribution by integrating it over the range of x. More importantly, dividing r/ ldln by 
P[r > ti\ yields the conditional density on the evidence x(ti) conditioned on no decision 
until time t\. 


4-1. Probability of hitting the lower threshold 

In trying to compute P[x(r) = —£], we view the two-stage process as two single-stage 
processes in sequence. Let T\ be the first passage time for the first stage by itself (without 
any deadline at t\) and define the random time f± = min{ri,H} which is a stopping time. 
Applying the optional sampling theorem to the exponential martingale for the first stage 
gives us 

exp(—2siXo) = E[exp(—2six(fi))] 

= E[exp(—2 six(t 1 ))|t 1 < ti]P[Ti < tf + E[exp(—2six(ti))|r 1 > H]P[n > H] 

= (exp(-2siC)P[x(Ti) = CCi < h] + exp(2siC)P[x(Ti) = —Cl^i < H])P[n < h] 

+ E[exp(—2six(ti))l(ri > tf)] 

= (exp(-2siC) + (exp(2siC) - exp(-2siC))P[a;(ri) = -C|n < H])P[n < h] 

+ E[exp(-2six(H))l(ri > H)], 


which yields 

„r / \ __> n . , n _ exp^SjSp) - E[exp(■ -2s 1 x(t 1 )) 1 {t x > Hj] - exp(-2siC)P[ri < H] 

Tl r i — i exp(2siC) — exp(—2siC) 

In this expression, both E[exp(—2s i a;(ti))l(r 1 > tf)] and P[ti < t-f may be obtained 
from ( [T7| ). 

The second stage is another single-stage process, this time starting at time t\ with a 
random initial condition x(t\) with distribution given by (17). Computing the expected 
value of the standard lower threshold hitting probability (|6j) with respect to the random 
initial condition Ad, i.e., x(ti) conditioned on T\ > t\ we obtain 

E[exp(— 2s 2 x(t 1 ))\Ti > ti] - E(-2s 2 C) 


[ z ( t ) = ~C\n > h] = 


exp(2s 2 C) - exp(-2s 2 C) 
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The E[exp(—2s 2 a;(ti))|r 1 > t^ term can be readily computed from Combining the 
previous two conditional expressions allows us to obtain 

P[x(r) = -C] = P[z(ti) = -c & n < ti] + P[x(r) = -C|ti > ti]P[ri > ti\. 

This probability depends on all of the parameters of the two-stage model: x 0 , pi, cp, /i 2 , cr 2 , (. 
4-2. Expected decision time 

To compute the expected decision time, we apply the optional sampling theorem to the 
martingale x(t) — pit with the stopping time f\ from the previous section to obtain 


xq = E[x(fi) - pi n] = E[x(ri) - HiT-^Tx < fi]P[ri < t x ] + E[x(ti) - pit^n > ti]P[ri > ti] 

= (C(l - 2 P[x(t 1 ) = -Cln < h]) - piEfrijr! < fi])P[ri < ti] + (E[:r(ti)|Ti > ti] - > h}. 

Solving for E[t1(ti < ti)] we obtain 

E [ rl ( ri < | a |j = (P[ r i < * 1 ] - 2P[x(r 1 ) = -C & n < ti])C - Xq + (E[a:(fi)l(ri > tQ] - /iAP[ri > fj) 

Much like in the previous section, we observe that the second stage is similar to a single-stage 
process starting at time t\ with a random initial condition determined by 0- Thus, the 
associated expected first passage time is 

Tn-r i + , (1-2 P[x(t) = -C|t >t 1 ])C-E[x(ti)|T 1 >ti] 

E[t|T! > t ij —ti-\ -. 

IE 

Combining the above expressions, we obtain 


E[r] = 


(P[ti < h\ - 2P[x(r 1 ) = -C & n < ti])C - x 0 + (E[a:(t 1 )l(r 1 > ti)] 


pi 


(1 - 2P[x(r) = -C| t > ti])P[ri > fi]C - E[x(t 1 )l(r 1 > ti)} 


IE 


4-3. First passage time density 

We now compute the first passage time probability density function. Let F(t; 6 , x 0 ) be 
the cumulative distribution function of the decision time for the single-stage process ([Tj) 
obtained by integrating (13). For t < ti, the two-stage model is identical to the first-stage 
model and the first passage time distribution is F(t] 9 1 , x 0 ). For t > ti, the first passage time 
distribution is 


F(t; 9 1 ,9 2 , x 0 ) = F(tp 6i,x 0 ) + E[F(t - tp 0 2 , x(ti))l(ri > ti)], 

i.e., the distribution function corresponds to trajectories that reach threshold before t\ and 
trajectories that reach threshold between t\ and t. The latter trajectories can be modeled as 
trajectories of a single stage process starting at time t\ with stochastic initial condition. The 
stochastic initial condition leads to the expectation operator on the second term. The dis¬ 
tribution of decision time conditioned on particular decisions can be computed analogously. 
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5. Analysis of the multistage model 


In this section we derive first passage time (FPT) properties of the multistage process 
defined in £2/2 using an approach similar to that employed throughout £j4j The model is 
viewed as n modified processes in sequence in which for each stage, the initial condition is a 
random variable and only the decisions made before a deadline are considered. For the i-th 
stage process with a known distribution of initial condition Xj_i, we derive properties of the 
FPT conditioned on a decision before the deadline C, along with the distribution of X t for 
i G {1,..., n}. This latter distribution, more precisely, is the distribution of x(fj) conditioned 
on the FPT for the i-th stage being greater than fj. We then use these properties sequentially 
for i e {1,..., n} to determine the FPT properties during each stage. Finally, we aggregate 
FPT properties at each stage to compute FPT properties for the whole multistage process. 
Our calculations have features similar to an idea of Diederich and Busemeyer (2006), who 
first proposed that the bias (i.e. initial condition) of a stage may have a time dimension. In 
some sense our formulas below elucidate how previous temporal stages of processing affect 
the bias of future stages. 


The extension of the FPT distribution computation from the two-stage model (Ratcliff 


1980) to the multistage model requires careful computation of expressions similar to (17) 


at the end of each stage. Also, in contrast to Ratcliff (1980), as in §J4|, our martingale 
based approach allows direct computation of probability of hitting a particular threshold 
and expected decision times. As in the previous section, this avoids integration of the first 
passage time density to compute these quantities. 

Throughout this section we use the following notations: 


• r = inf{f > 0|x(t) ^ (—CC)}j the first passage time through either threshold for the 
entire multistage process; 

• Tj = t\t > ti- 1 , the first passage time for the f-th stage ([3]) without any deadline; 

• 0i = and 0 1:e = (/zi,..., oi,..., erg, ti ,..., tg, () representing the pa¬ 

rameters for the i-th stage and stages 1,... ,£, respectively; 

• Si = /Zj/of, the z-th stage ratio of signal to squared noise. 


Here we are concerned with computation of FPT properties and allow 0\- n to be free 
parameters. For scenarios such as estimation of parameters, to ensure identifiability of the 
parameters, all diffusion rates may be set equal to unity. However, such cases are beyond 
the scope of this manuscript and we do not discuss these issues here. 


5.1. FPT properties of the i-th stage 

For the i-th stage, the initial condition A7 i _ 1 is a random variable and only decisions 
made before the deadline t. t are relevant. The analysis thus focuses on the random variable 
Aj_i and the random time T t . Conditioned on a realization of AWi, the density of A 7 * can 
be computed using ( Jl7| ). If the density of Aj_i is known, then the unconditional density 
of Xi can be obtained by computing the expected value of the conditional density of X t 
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with respect to X t _ L . Since the density of X 0 is known, this procedure can be recursively 
applied to obtain densities of for each i e {1, Formally, the joint density 

gf dln (x,r]Xo,9i : i) of the evidence x(U) and the event t >U is 


gf dln (x, n ; x 0 , U, 0 1:i ) = -^-P [x(ti) < x & n > t { \ = [g ddln (x, n - Vi; - t*_i, 9i)}, 

ax 

(18) 


where denotes the expected value with respect to X t _\. Note that 


P [x(ti) < X & Ti > ti\ = P [x(ti) < X 


Ti > ti]P[Ti > ti] = P [Xi < x]P[Ti > ti). 


Thus, the density of X i: i.e., x{ti) conditioned on r t > U is determined by dividing gf dln by 
P[rj > ti\ which can be computed by integrating g ddln over the range of x(f*). Note that 
the parameters in g 4 ddln are Xo and 9 \ : jj this highlights the fact that the distribution of X t 
depends on all previous stages. 

Similarly, the FPT density for the i-th stage conditioned on a realization of X t _\ can 
be computed using (13), and the unconditional density can be obtained by computing the 
expected value of the conditional density with respect to AVn 


fi(t] Xq, 0 1:i ) = ^P [Ti <t]= E AVl [/(t - ti_ i; Xi_ 1 , Oi)], (19) 

where t > £j_i. The cumulative distribution function F r (t: xq, 9 i :i ) = P[rj < t] is obtained by 
integrating fi(t; xq, 9 i ;i ). Note that every trajectory crossing the decision threshold before U 
does so irrespective of the deadline at ti. Thus, the expression for density does not depend 
on ti . 

To conclude this section, we state equations for computing hitting times, mean decision 
times, and first passage time densities, conditional on a response during the i-th stage. The 
derivation of these expressions are found in |Appcndix B 

(i) The probability of hitting the lower threshold given that a response is made during the 
i-th stage, denoted by P^(x 0 ,9 1:i ) := 1 - T) + (x 0 , 9 1:i ) := P[x(r) = ~C\U-i < r < U], is 
given by 


P~ 


' Ex I _ 1 [exp(—2s i .Y i _i)] - E Ai [exp(-2s i jy i )]P[r i > U\ - exp(-2siC)P[n < U] 
(exp(2siC) - exp(-2siC))P[n < U\ 

1 (E AVl [AVi] - E. Vi [A]Ph > U]) 

. 2 2CP[n < U] 


if fit ± 0, 

( 20 ) 

if fii = 0. 


P+(x 0 , 9i : i) is the probability of hitting the upper threshold during the f-th stage. These 
expression depends on all of the model parameters up to and including the i-th stage 
(i.e. Pp := P~ (x 0 , 9i :i ) and P+ := Pj + (x 0 , 0 1:i ) both depend on x 0 and 9 i ;i ) and will 
be used in subsequent calculations. 
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(ii) The joint FPT density for the i-th stage process and a given upper/lower response, 
denoted by jf(U x 0 , 6 l:i ) = )|P[r* <tk x(r*) = ±C], is given by 

fi{t\x o,0 1:i ) = [Ti <tkx{Ti) = C] =P Xi _ 1 [f + (t (21) 

f~(t-,x 0 ,e 1:i ) = ^P [Ti < tkx(Ti) = -C] = -U-i; Xi-I,9i)], ( 22 ) 


(hi) 


where the functions f ± (t ; xq, Oi) are taken from (14) and 0 - Again, these expressions 
depend on all of the multistage model parameters up to and including the i-th stage 
(i.e. both fi(t;xo,9i:i) depend on xq and 9\ : i). 

The mean decision time given a response during stage i, denoted by mDTj(xo, 0i : *), is 
given by 


mDTj(x 0 , 9\ : i) = E[tj|T j < U] 

(1 - 2 Pr)(P[Ti < ti] - E^.JAVi] + E x .[Xj]P[Ti > U] - - ij_i)P [n > tj 


= 


ti-1 + 


ti-1 + 


T, < ti] 

C 2 P[n < ti] - + E^[A7]P [ Ti > u] - a\{ti - h-!)P [Ti > ti] 

of Phi < ti] 


if ^ ^ 0, 

if ^ = 0. 

(23) 


lv 


The mean decision time conditioned on a given upper/lower response made during the 
i-th stage, denoted by mDTf (x 0 , 9 1:i ), is given by 

mDT+(z o ,0i:i) ^P\Ti\x{ji) = C & Ti < ti] 
mDT 4 + (x 0 ,9 1:i ) 


P+(x 0 ,9 1:i )P[Ti < ti] 

1 


P+(xo,9i : i)P[Ti < ti] 


ti-iP[x(Ti) = C] +Exi_! mDT (A"j_i, 9i) 


- E 




mDT (Xi, 9i ) - tiP[x(ji) = C |Ti > U] IP [n > U] 


(24) 


mDT; (xo,0i ;i ) =E[r i |a;(r i ) = ~C & A < U] 
mDT, : (x 0 ,9 l:i ) 


Pi (Xq, 01:i)P[Tj < ti] 

1 


Pi {XQ, 9i : i)P[Ti < ti] 


ii_iP[x(rj) = -C] + E Xi _! mDT (A,_i,0*) 


E 


A; 


mDT (X if 9i ) - tiP[x(Ti) = -C| Ti > ti]jP[Ti > ti]j , 

(25) 
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where mDT i (x 0 , 0 1:i ) = E[r,;l(x(rj) = ±£ & t* < L)] and is calculated using (11) and (12), 
P[z(Ti) = ±C] =E Xi _ 1 [P ± (X i . 1 ,e i )\, and Pfcfa) = ±C | A > U] = E Xi [P±pQ, Of)] is 
calculated using (|6|) and (18). 


5.2. FPT properties of the multistage model 

For a given multistage process (J2| with initial condition xq, we sequentially compute all 
of the distributions of the initial conditions X t for each i G {1 using (18). Then, 

we compute the properties of the FPT associated with the i-th stage. Finally, the total 
probability formula aggregates these results into FPT properties of the entire multistage 
model. The calculations are contained in Appendix C and are given in terms of earlier 
formulas. In the following, we omit the arguments of functions whenever it is clear from the 
context. 

(i) Let t > 0 be given such that t e (i fc _i, f fe ] for some k e {1,..., n}. The FPT distribution 
for the multistage process is 


k—1 k —1 

P[t < t] = 1 - P [ri > ti] + P[r fc < t] P [Ti > ti\. (26) 

i— 1 2=1 

Note that n*=i E[p > U] — Ppr > tk- i] an d P[T fc < t] — P[r < t\r > 

(ii) The mean decision time, denoted by mDT ms (x 0 , 9i :n ), for the multistage process is 

n i— 1 

mDT ms (x 0 ,9i :n ) = E[r] = ^ (E[Ti|7i < 2*]P[Tj < U] J^Pfa > tf}) 

i=l j =1 

n z—1 

= ^2 Oi :i )P[Ti < U] P [tj > U}'). (27) 

i=l j =1 

Put simply, the expected decision time is the sum of the expected decision times for 
the individual stages (mDTj) weighted by the probability of the decision in each stage 

(p[tj < ti\ n;;! p[p >u]). 

(iii) The probability of hitting the lower threshold , denoted by P~ s (x 0 ,9i. n ), is 


2 — 1 


P ms(XO, dim) = 1 - ^s(A, 0 1:n ) = ^ (P, P fa < U] JJPfo > tj 

3 = 1 


(28) 


2=1 


This expression is similar to (27), with the probability of hitting the lower threshold be¬ 
ing the sum of the hitting probabilities for each stage (P“) weighted by the probability 
of the decision in each stage (P[r,; < u] n;J, pfo > u]). 
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(iv) The mean decision time conditioned on hitting the upper/lower threshold is 


Z—1 


E[r|x(r) = C] = -pjr E (E[ri|x(ri) = C & n < U] P+ P[p < U] JJPfo > U]j (29) 


3 =1 
l—l 


i=1 


E[r|x(r) = -C] = -pr E [ E [ T i\ x ( r i) = C &Ti<ti] P i P [n < ti] JJPfa > U]). (30) 


3 = 1 


Note that E[rj|a:(rj) = C & T i < U]Pi = E[rjl(x(rj) = £)|ri < ti] and P[rj < 


ti] nt=iP[ r j > U] = P[£j-i < r < ti]. Similar to (27), these equations show that 


the conditional decision times are weighted sums of the expected conditional decision 
times for each stage, with the weights being the conditional probability of the decision 
in each stage. 

(v) The FPT cumulative distribution functions conditioned on hitting upper/lower thresh¬ 
old are 


1 


k -1 


[r < t\x(r) = (] = — I P[r fc < t & x(r k ) = C] JJ P [ r i > 

■*ms \ j = i 

k —1 z—1 


E p [n < u & x{ji) =c] XT > ^ 


i=l 


3 =1 


(31) 


1 / fc_1 

[r < t\x(r) = -C] = — ( P(r fc < t & x(r fc ) = ~C) JJPfo > A?] 

ms \ j = i 


k— 1 z—1 

E P [ T * - ^ & X ( X i) = _ C] > tj 

3 =1 


z=l 


(32) 


Note that 11^=1 Pfa > £?] = P[t > £*-i] and P[r* < ti & x(rj) = C]P[r > ^-i] — 
P[L-i < r < Zj & x(r) = (]. 


6. Time-varying thresholds for the multistage process 

The results in §|5]were obtained under the assumption that the thresholds are constant 
throughout each stage. Now suppose that the thresholds for the i-th stage are ±Q, he., 
piecewise constant thresholds. If the upper thresholds decrease at time t t (i.e. C*+i < Ci) 
and x{ti) is in the interval (£j + i,£j), then the path is absorbed by the upper boundary, and 


the probability of this instantaneous absorption is calculated by integrating (18) from £ i+ i 
to (i . Likewise, the probability of instantaneous absorption into the lower threshold at ti 


is determined by integrating (18) from — Q to —0+i- The density of X t is then found by 


truncating the support of the density in (18) to (— Ci+i, C*+i) an d normalizing the truncated 


density. In the cases where the upper threshold in the (i + l)-th stage is larger than the 
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Figure 1: Illustration of the key ideas for computation of FPT properties of the multistage model with 
piecewise constant drift and thresholds. The distribution of evidence x{ti) conditioned on no decision until 
U serves as the distribution of initial condition for the *-th stage process. If threshold £i+i < Q, then 
probability of instantaneous decision at time ti is computed as the probability of x[ti) conditioned on no 
decision until ti not belonging to the set (—If threshold (j + i > Q, then there is no instantaneous 
decision and only the support of x{ti) conditioned on no decision until ti is increased to (— Ci+i> C*+i)- 


upper threshold in the i-th stage, i.e., Ci+i > G, there is no instantaneous absorption, and 
the density of X t is found by extending the density in (18), assigning zero density to the 
previously undefined support (see Figure [l|. In all cases, the new, updated X t may be 
used for computations dealing with the (i + l)-th stage of the multistage model. Codes 
implementing all of the formulas through §|5] with time-varying thresholds may be found at 


https://github.com/sffeng/multistage. In Appendix D ( we describe how these ideas 


extend to a time varying Ornstein-Uhlcnbeck (0-U) model. 


7. Numerical examples 

In this section we apply our calculations from f|5] and '|6] to a variety of numerical ex¬ 
periments. In doing so, we compare the theoretical predictions obtained from the analysis 
in this paper with the numerical values obtained through Monte-Carlo simulations, thereby 
numerically verifying our derivations above. We also provide examples illustrating time pres¬ 
sure or changes in attention over the course of a decision process, and demonstrate how our 
work can help to find the optimal speed-accuracy trade-off by maximizing reward rate (or 
any other function of mean first passage time, threshold-hitting probability, and reward). 
Unless otherwise noted, Monte Carlo simulations were obtained using 1000 runs; relatively 
few runs are used so that curves are visually distinguishable. Stochastic simulations use the 
Euler-Maruyama method with time step size 10 -3 . All of the above calculations have been 
implemented in MATLAB, and all codes used to produce the figures in this section may be 
found at https://github.com/sffeng/multistage. 
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7.1. A Four-stage process 

Consider a four stage process with drift and diffusion rates given by (//i, p 2 , 1 ^ 3 , Ta) = 
(0.1,0.2,0.05,0.3) and (<7i, <t 2 , cr 3 , cr 4 ) = (1,1.5,1.25,2), respectively, with (t 0 ,ti,t 2 ,t 3 ) = 
(0,1, 2, 3) and initial condition xq = —0.2. The cumulative distribution function (CDF) 
of the unconditional and conditional decision time for ( = 2 obtained using the above 
analytic expressions (solid lines) and Monte-Carlo simulations (dotted lines) is shown in 
Figure 2(a) Similarly, the unconditional mean decision time, the lower hitting probability, 


and the mean conditional decision times (for upper/lower responses) are shown in Figure [2(b)| 
as a function of threshold (. Note that the analytic expressions match closely with quantities 
computed using Monte-Carlo simulations. Also, notice that the CDF almost looks like a 
double-sigmoidal function (it starts to saturate around 0.6 before picking up and eventually 
saturating at 1) due to the drop in drift rate from 0.2 to 0.05. 

7.2. FPT Distribution for a process with Alternating Drift 

I 11 this example, the sign of the drift rate changes from stage to stage. This may be 
used to describe situations in which evidence accumulation changes dynamically with the 
decision-maker’s focus of attention. For instance, Krajbich et al. (2010) have shown that 


the process of weighing two value-based options (e.g., foods) can be modeled with a process 
in which drift rates vary based on the option being attended at any given moment. We 
consider such a case using a 30-stage model in which the drift rates 1 and —0.75 alternate 
(i.e., fii — 1, H 2 — —0.75, /r 3 = 1, ...) to capture a situation in which the decision maker’s 
attention alternates between two options, one of which has greater perceived value (higher 
drift rate) than the other. Let to = 0 and the remaining 29 stage initiation times be a fixed 
realization of 29 uniformly sampled points between 0 and 10. Assume xo = 0, ( = 2, and 
let the diffusion rate be stationary and equal to unity (<7; = 1). The unconditional and 
conditional FPT distributions in this scenario obtained using both the analytic expressions 
(solid lines) and Monte-Carlo simulations (dotted lines) are shown in Figure [3j Note that the 
analytic expressions match closely with quantities computed using Monte-Carlo simulations. 


7.3. FPT Distribution for a model with gradually time-varying drift 

Changes in evidence accumulation may occur gradually over time. For instance, White 
et al. (2011) proposed a “shrinking spotlight” model of the Eriksen Flanker Task, a task 


in which participants responding to the direction of a central arrow are influenced by the 


direction of arrows in the periphery (see also Servan-Schreiber et al., 1990 Liu et al., 2009 


Servant et al. 2015). According to these models, evidence accumulation is initially influenced 


by all of the arrows (central as well as flankers, which may drive an incorrect response, 
modeled here as as a lower threshold response) but as the attentional spotlight narrows the 
drift rate is gradually more influenced by the central arrow alone. The multistage model, 
in spite of having discontinuous changes in parameters, can still be used to approximate a 
model with gradually time-varying parameters. 

As a demonstration, we use a 20-stage model as an approximation to a model with con¬ 
tinuously time-varying drift rate. Assume Uj = 1, xq — 0, £ = 2, and let the stage initiation 


19 
























(a) Unconditional and Conditional First Passage Time Distribution 



(b) Unconditional Mean Decision Time, Lower Hitting Probability, and Con¬ 
ditional Mean Decision Time 

Figure 2: FPT calculations for a four-stage process with drift rates (/zi, fj, 2 , M 3 ,/h) = (0.1,0.2,0.05,0.3), 
diffusion rates ( 01 , < 72 , cr 3 , 0 - 4 ) = (1,1.5,1.25,2), stage initiation times {to,ti,t 2 ,t 3 ) = (0,1,2,3), and initial 
condition xq = —0.2. The FPT distribution is computed for threshold ( = 2. 
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Figure 3: Unconditional and conditional FPT distributions for a 30-stage model with alternating drift rate. 
The drift rates are (/zi, /r 2 , /r 3 ,...) = (1, — 0.75,1,...), diffusion rate at each stage is unity, the threshold 
( = 2, and stage initiation times are equally spaced throughout the interval (0,10). 


times t 0 ,ti,... ,t 19 be equally spaced throughout the interval [0,5]. Furthermore, suppose 
the drift rate during the i-th stage is —0.2 + 0.0263(i — 1). The unconditional and conditional 
FPT distributions for such a 20-stage process obtained using the analytic expressions (solid 
lines) and using Monte-Carlo simulations (dotted lines) are shown in Figure [fj Note that 
the analytic expressions match closely with quantities computed using Monte-Carlo simula¬ 
tions. We also show the error due to piecewise constant approximation of the drift rate as a 
function of number of stages (Figure [dj right), ft should be noted that even for 5 stages the 
approximation error is very small. 


7-4- Collapsing Thresholds 

One may be interested in modeling a decision process in which thresholds are dynamic 
rather than constant across stages. This can be used to describe discrete changes in choice 
strategy, or a continuous change in thresholds over time; the latter approach has been success¬ 
ful at describing behavior under conditions that either involve an explicit response deadline 
(e.g., Milosavljevic et al. 2010; Frazier and Yu[ 2008) or where there is an implicit opportu¬ 
nity cost for longer time spent accumulating evidence (Drugowitsch et al., 2012). Recently, 


(Voskuilen et ah, 2016) used analytic methods to model collapsing boundaries in order to 


compare fixed boundaries against collapsing boundaries in diffusion models. 

Here, we model such a situation using a 20-stage process, as an approximation to a 
diffusion model with continuously collapsing thresholds ,i.e., ( 0 with time, the drift rate 

and the diffusion rate are constant and equal to 0.15 and 1, respectively, x 0 = 0, and stage 
initiation times t 0 , • • •, tig are equally spaced throughout the interval [0, 5]. The threshold 

in the i-th stage is ^ = 3 —^(i —1). The unconditional and conditional FPT distributions for 
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Figure 4: Unconditional (left) and conditional (middle) FPT distributions for a 20-stage process with grad¬ 
ually increasing drift rate. The drift rate for the i-th stage is /u, = —0.2 + 0.0263(i — 1), diffusion rate at each 
stage is unity, the threshold ( = 2, and stage initiation times are equally spaced throughout the interval [0, 5]. 
Right: The total squared error between the simulated and analytic CDFs (left) decreases with increasing 
discrete stages of the model. 10,000 simulations were used for each approximation. 


such a 20-stage process obtained using the analytic expressions (solid lines) and using Monte- 
Carlo simulations (dotted lines) are shown in Figure [5] Note that the analytic expressions 
match closely with quantities computed using Monte-Carlo simulations. 


7.5. Optimizing the Speed-Accuracy Trade-off in a Two-stage Model 

Human decision-making in many two alternative forced choice signal detection tasks has 
been successfully captured by the single stage model. In such tasks, hitting the upper/lower 
boundary is interpreted as a correct/incorrect response. The accuracy of a decision can 
then be determined by the sign of the drift rate - if /i is positive, participants are said 
to be more accurate the more likely they are to hit the upper threshold and more error- 
prone the more likely they are to hit the lower threshold. One then assumes without loss of 
generality that the drift rate is positive, in which case the lower hitting probability is called 
the error rate and the upper hitting probability is called the accuracy. Here, the choice 
of threshold dictates the speed-accuracy trade-off, i.e., the trade-off between a fast decision 
and an accurate decision. In the previous examples, the thresholds have been known and we 
have characterized the associated error rate and first passage time properties. These can be 
used to define a joint function of speed and accuracy that may dictate how humans/animals 
choose to set and adjust their threshold. In particular, it has been proposed (Bogacz et ah 


2006) that human subjects choose a threshold that maximizes reward rate (RR), defined as 


RR = 


1 -P~ 

E[r] + T n( j 


(33) 


where T n( j is the sensory and motor processing time (non-decision time) and P~ and E[r] 
are computed using the expressions derived in §|5j 

The reward rate for a two-stage process is shown in Figure [6} For the set of parameters 
in Figure [6j reward rate is maximal at approximately (Ci 5 C 2 ) = (0.06,0.01). Thus, the 
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Figure 5: Unconditional and conditional FPT distributions for a 20-stage process with collapsing thresholds. 
The drift rate and diffusion rate at each stage are 0.15 and 1, respectively. The stage initiation times are 
equally spaced throughout the interval [0, 5], and the threshold at i-th stage is Q = 3—^, for * = {1,..., 19}. 


maximizing reward rate in this setting interestingly requires that the threshold across stages 
be different. 

Setting the threshold to be constant across stages (£i = £ 2 ), we can compare how reward 
rate changes with this constant threshold in a single-stage (traditional) versus a two-stage 
mode. As shown in Figure 7(a), we find that this reward rate function is unimodal in a single- 


stage model (as previously observed) whereas it is bimodal in a two-stage model. Figure |7(bJ 
explores this parameter space in greater depth and shows that the curvature of reward rate 
(and in particular the relative height of its first and second modes) vary as a function of 
the length and drift rate of the first stage (for example). As a result, this analysis reveals a 
discontinuous jump in optimal threshold as these parameters vary. Whether individuals are 
sensitive to these discontinuities when setting thresholds for a multistage decision-making 
task deserves further exploration. 


8. Discussion 

In this paper we analyze the first passage time properties of a Wiener process between 
two absorbing boundaries with piecewise constant (time-dependent) parameters, which we 
call a multistage model or multistage process. Our main theoretical results, collected in §[5j 
add to previous work on analyzing time-dependent random walk models in psychology and 
neuroscience. Broadly speaking, these can be split into three approaches. One approach is 


the integral equation approach introduced and developed in Smith (1995, 2000); Smith and 


Ratcliff (2009). Another approach is the matrix based Markov Chain approximation which 
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Figure 6: Reward rate as a function of the thresholds for the two-stage model. The parameters are pi = 0.5, 
H 2 = —0.1, ai = cr 2 = 0.1, x 0 = Ci/2, and t\ = 0.1. The non-decision time T nd = 0.3. Note that the 
maximum reward rate is achieved for ^ C,-2- 
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(a) Reward rate for single- and two-stage models 


(b) Reward rate for two-stage model 


Figure 7: (a) For the single-stage model H\ = 0.2, o\ = 0.1 and Xq = 0, while for the two-stage model 

/Lti = 0.1, /x 2 = 0.5, or = (72 = 0.1, Xo = 0 and t\ = 0.15. The non-decision time T nd = 0.3. The reward rate 
for the single-stage model is a unimodal function and achieves a unique local maximum, while the reward 
rate for the two-stage model has two local maxima, (b) Optimal threshold for two-stage model obtained 
by maximizing reward rate. The left panel shows the variation of the optimal threshold as a function of t\ 
and /i-|. The other parameters are fi 2 = 0.5, <J\ = <12 = 0.1, and x 0 = 0. The regions of the contour plot 
associated with t\ = 0 and fi 1 = 0.5 correspond to the single-stage model. 


has been applied to a wide variety of multi-attribute choice settings (Diederich and Buse- 


meyer[ 2003 Diederich and Oswald 2016 2014). A third approach analyzes the backward 


partial differential equation associated with the multistage process (Ratcliff, 1980 Heath 
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1992). The results of fjs] align most closely with the work of the third approach, as we also 


directly analyze the multistage stochastic process, albeit with different techniques. Whereas 
Ratcliff (1980) and Heath ( [199 2) analyzed a multistage process by solving the Kolmogorov 
partial differential equation, we employ martingale theory (e.g. the optional sampling the¬ 
orem) in order to obtain analytic results which extend those of previous studies. In doing 
so, our work also builds on martingale-based analyses described by Smith (|1990) for a single 
stage model. 

The modeler utilizing time-dependent random walk models in decision making should 
be aware of all of the above approaches, as one approach may demand additional approxi¬ 
mations compared to another, due to differences in how they discretize temporal dynamics. 
For example, when modeling experiments with continuously (gradually) changing stimuli 
(e.g., White et al. (2011); Zhang et al. (2014)) one may find the integral equation approach 
more natural, given the continuity and smoothness assumptions built in to the techniques. 
Diederich and Oswald’s model also allows for continuous changes in boundaries. In contrast, 
our methods must approximate the underlying gradually changing drift with a piecewise 
constant function, thereby inducing some additional error in the calculations of first passage 
times. If, however, the application lends itself well to discrete changes in drift rate or thresh¬ 
old (e.g., Krajbich et al. (2010)), then both the matrix based methods and methods discussed 
in this paper may be more natural, as they explicitly consider such discrete changes in the 
underlying calculations. Ornstein-Uhlcnbeck processes lend themselves more naturally to 
the matrix based approach of Diederich and Busemeyer (2003), compared to our analysis in 


Appendix D, since our analysis requires a change in time scale that the matrix method does 


not. In the specific case of a multistage model ([2]), our work provides semi-analytic formulas 
that can be easily computed or studied further (see below). In general, the modeler should 
be ready to employ the most suitable approach given their situation. 

The reader should also be aware of the software tools available for each approach, as 
there already exist several good non-mart ingale based software packages for computing FPT 
statistics. One package implementing the integral equation approach is that of Drugowitsch 


(2014), which computes first passage time densities using the stable numerical approxima¬ 
tions developed in Smith (2000). More recently, highly optimized codes for a broad class 
of diffusion models have been developed by Verdonck et al. (2015), with implementations 
on both CPUs and GPUs. Compared to other available codes, Diederich and Busemeyer 


(2003)’s matrix approach is substantially simpler and more elegant to implement - one can 
compute desired choice probabilities and mean decision times in less than a few dozen lines 
of MATLAB code. For practitioners wishing to write all of their own model code from 
scratch, this may be a considerable advantage. The MATLAB code released with this re- 
portj^] provides implementations of the results from |5] and allows one to reproduce all of 


the figures from §[7| Unlike the work of Verdonck et al. (2015), our work is not immedi¬ 


ately focused on developing a rapid numerical tool for simulation, but rather on introducing 
martingale theory as a useful mathematical tool for analyzing multistage decision models. 


’https://github.com/sffeng/multistage 
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Thus, the codes released with this report are not intended to compete with the efficiency of 
the aforementioned codes, which have been highly optimized and timed for throughput, but 
instead to demonstrate the simplicity and effectiveness of our analysis. These considerations 
notwithstanding, our results do suggest promising avenues for future numerical work. Par¬ 


ticularly relevant is work by Navarro and Fuss (2009); Blurton et al. (2012); Gondan et al. 


(2014) who developed efficient numerical schemes for evaluating the relevant infinite sums 
involved in FPT calculations. Similar methods could be applied to results in 15.1 and 15.2 to 
develop efficient multistage codes, which could in turn contribute to the growing collection of 
numerical tools available for practitioners using diffusion models to study decision making. 

In 1(3] we re-derived classical mean decision times and choice probabilities equivalent to 


results first derived for a discrete time random walk model using the Wald identity (Laming 


1968 

Link and Heath, 

1975 


Link, 

1975 

Smith 

1990 

), which itself is a corollary of the 

optional sampling theorem. 

Smith 

(1990 

i notes that results concerning the discrete time 


random walk model with Gaussian increments and the continuous time single stage Wiener 
diffusion model of §|3] should be equivalent, which is indeed the case. Furthermore, the 
moment generating function derived in 13.3 is identical to that obtained by (Smith 1990, p. 


9) via the Wald identity. Our aim in presenting 1|3] was to introduce the martingale analysis 
used in l|5] by first demonstrating its utility in re-deriving these classical results. We hope 
these calculations provide an intuition for the computations in l|5] in a less technical setting. 
For reviews of the classical results on random walks in psychological decision making, see 

328-334). 


(Townsend and Ashby, 1983 

pp. 300-301) and ( 

Luce, 

1986 


pp. 


One particularly useful aspect of equations (20) through (32) is that they demonstrate 
how various other measures of performance, such as the lower-threshold hitting probability 
during each stage, evolve as the underlying dynamics change. Using these, one may effi¬ 
ciently compute a variety of behavioral measures of performance without resorting to first 
computing the FPT densities. Our results may also serve as a starting point for further anal¬ 
ysis of more complicated stochastic decision models. In appendix D, we show how our results 
apply to Ornstein-Uhlcnbeck processes, which approximate leaky integration over the course 


of evidence accumulation, e.g., the Leaky Competing Accumulator model [LCA] (Usher and 


McClelland, 2001). Given that the LCA itself can in certain cases approximate a reduced 


form of more complex and biologically plausible models of interactions across neuronal popu¬ 
lations (e.g., Wang, 2002 Wong and Wang, 2006 Bogacz, 2007), our work may help analysts 
better understand time-varying dynamics within and across neural networks, and how such 
dynamics relate to complex cognitive phenomena. 

Beyond its analytic and numerical utility, our uses of martingale theory, and the optional 
sampling theorem in particular, provide some theoretical insights into random walk decision 
models. For example, a single stage model with zero drift rate (i.e. Wiener process) and 
lower/upper thresholds at 0 and 1 necessarily has the qualitative property that the proba¬ 
bility of hitting the upper threshold equals the initial (nonrandom) position xq . Although 
these results are known in the probability literature (Doob, 1953), they provide practitioners 
and experimentalists a unique way of planning and analyzing experiments based on decisions 
thought to evolve according to diffusion processes. In summary, we hope that the tools de- 
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scribed in this article and the results of 1|5] will encourage computational and mathematical 
analyses of decision models involving time-dependent parameters. 
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Appendix A. Alternative expression for FPT density of the single stage model 


An alternative expression to the FPT density (13) that can be obtained by solving the 
Fokker-Planck equation (Feller, 1968) is: 
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Appendix B. Derivation of expressions in §5.1 


We first establish (20). First consider the case ji l > 0. Let be the filtration 

defined by the evolution of the multistage process ([2]) until time t conditioned on r > t^i. 
A filtration can be thought of as an increasing sequence of available information. For some 
s E (L-i, t), as shown in fj3j {exp(— 2s i x(t))} t > t ._ 1 is a martingale, i.e., E[exp(— 2s i x{t))\J ri f \ = 
exp(— 2six(s)). 

Furthermore, fj := min{rj,tj} is a stopping time. Therefore, it follows from optional 
sampling theorem that 


E [exp (—2s j Aj_ i)] = E[exp(-2six(fi))] 

= E[exp(-2s i a;(r i ))|r i < f*]P [n < U] + E[exp(-2s i A 7 'j)]P[r i > ti\] 

= (exp(-2siC)(l - P~) + exp( 2 sj£)Pj“)P[Tj < U] + E[exp(-2sjAj)M T * > ^]]- 

Solving the above equation for P~ , we obtain the desired expression. 

For /x,; = 0, , as shown in §[3j {x(t)} t > ti is a martingale. Therefore, applying the optional 
sampling theorem, we obtain 


E[A;_i] = E[^(a)] 

= E[x(ri)|ri < ti]P[ri < U\ + E[A"j]P[xj > U}\ 
= (1 - 2P~)(P[Ti < ti\ + E[^Q]P[p > ti\. 


Solving the above equation for P ~, we obtain the desired expression. 

The formulas (21) and (22) immediately follow from applying expectation to (14) and (15), 
respectively. 


To establish (23) for //,- ^ 0, we note from §|3]that for the i-th stage, {x(t) — is a 

martingale. Therefore, applying the optional sampling theorem, we obtain 


E[Aj_i] - fiiti _i = E[x(fj) - HiTi\ 

= E [x(Ti) - iliTi\Ti < ti]P[Ti < ti] + E[(Aj - liiti}P[Ti > ti] 

= (C(l - p ~) - c p i - < L])P[p < ti] + (E[Aj] - fj,iti)F[Ti > ti}. 

Solving the above equation for E [r t |r,; < ti] yields the desired expression. 

For Hi = 0, we note from f|3]that {x(t ) 2 — l is a martingale. Therefore, applying 

the optional sampling theorem, we obtain 

EM - <7^, = E[x(f) 2 - a 2 fi] 

= E [x(Ti) 2 - a 2 Ti\Ti < L;]P [Ti < U] + E[A 2 - a 2 ti]F[Ti > U]] 

= (C 2 - cr 2 E[Ti\Ti < ii])P[Ti < ti] + (E[A 2 ] - (T 2 ti)P[Ti > ti}. 

Solving the above equation for E[r, : | r t < ti] yields the desired expression. 
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Next, we need to establish that the Laplace transform of the density for the FPT for a 
particular decision made before t t is 


E[exp(-ar i )l(x(r i ) = Cl n < U)} 


exp(-atj_i)E[TC lgf (/i», cn , C, Aj-i, op] ~ exp(-a^)E[T° lgl (^, a h C> x h oQMt > U}\ 

P [n < U] 


mgf / 


E[exp(-ar i )l(a;(r i ) = —C Vi < U)] 

_ exp(-atj_i)E[r mgf (/ii, <7j, C, Xj_i, a)] - exp(-o^)E[r mgf (/q, cq, C, X u aQPfc > fj] 

Pfo < U] 

To establish this, we consider the stochastic process exp(Ax(t) — A/Xj t — X 2 a 2 t/2). From §[3j 
we note that (exp(Ax(f) — A fj.it — X 2 a 2 t/2)} t > til is a martingale for each A G M, i.e., 

E[exp(Ax(t) — A [lit — \ 2 a 2 t / 2)\Fl] = exp(Ax(s) — A/i^s — A 2 cx 2 s/2). 

We choose two particular values of A: 

A = ~Vi ~ y/Vi+ 2aa i and A _ ~Vi + y/rf + ZaVi 
a 2 ’ an 2 a 2 

Note that for A G { Ai, A 2 }, A/ijf+A 2 cr 2 f/2 = a. Therefore, stochastic processes (exp(Ai:r(f) — at)} t > 0 
and (exp(A 2 a;(t) — at)} t > ti l are martingales. Now applying the optional sampling theorem, 
we obtain 


E[exp(AiW-i - oiti_ 1 )] = E[exp(Aix(f) - erf)] = exp(AiC)E[exp(-ari)l(a;(ri) = C & n < U )] 
+ exp(-AiC)E[exp(-ar)l(x(ri) = ~C & n < U)} + exp(-at i )E[exp(AiW)]P[^ > U\. 

(B.l) 

Similarly, 

E[exp(A 2 X;_i - oiti- 1 )] = exp(A 2 C)E[exp(—arj)l(a;(r) = C & P < U)\ 

+ exp(-A 2 C)E[exp(-arj)l(a;(r) = ~C & n < U)} + exp(-at i )E[exp(A 2 X i )]P[r i > £*]]. 

(B.2) 


Equations (B.l) and (B.2) are two simultaneous equations in two unknowns E[exp(—ar*) 1 (x(r,j) = 
C & r,; < tj)] and E[exp(— aTj)l(x(ri) = —( & t* < tj)]. Solving for these unknowns, we 
obtain 


E[exp(-ar i )l(a;(r i ) = C & n < U )] = 


exp(-at i _ 1 )E[exp(Ai(AVi + 0) 


exp(2AiC) - exp(2A 2 C) 

- exp(A 2 (Aj_i + C))] - exp(-at ; :)E[exp(Ai(Ai + ()) - exp(A 2 (Aj + C))]P[p > U] ), 
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and 


E[exp( —Q!Tj)l(>(Ti) = -C & n < ti)} = 


exp(-atj_i)E[exp(-Ai(C - ^Q-0) 


exp(-2AiC) - exp(—2A 2 C) 

— exp(—A 2 (C - AA-i))] -exp(—ati)E[(exp(-Ai(C - JQ)) -exp(-A 2 (C - A' i )))]P[r i >U]). 

Simplifying these expressions, we obtain the desired expression. 


Finally, (24) and (25) follow from differentiating the Laplace transform with respect to 


—a, and then evaluating at a = 0. 

Appendix C. Performance metrics for the overall multistage process 


We start by establishing (26). Since t e (4-i,4]> 

P[t < t] — P[r < t & r < 4-i] + P[r < t & r > 4-i] 
= P[r < 4- i] + P[t < t\r > 4_i]P[r > 4- i] 


fc-i 


fc-i 


1 _ II P t r > ^l r > + F lw < t] JJP[r > 4|r > t 

2=1 2=1 
fc— 1 fc—1 

1 - n p [ r * > + p [ Tfc - ^ n p [ r * > 


*—ij 


2=1 


2=1 


We now establish (27). We note that 


E[t] = El E[rl(4-i < r < ti)] 

2=1 

n 

= y^E[rl(r < ti)\r > 4-i] F b" > 4_i] 


2=1 


2—1 


2=1 


E ( E N r i < ti]^[Ti < ti] J^P[Tj > tj]j. 

3 =1 


To establish (28), we note that 

n+1 

P~ = P(x(r) = —C and 4-i < r < ti 


2=1 

71+1 


E] P(x(t) — C and r < U\r > 4-i)P[t > 4-i] 


2=1 

n 


2 — 1 


^(i’-p[r i <i i ]nPfe>i j ]). 

3 = 1 


j=l 


similarly to (26). 


Equations (29) and (30) follow similarly to (27), and Equations (31) and (32) follow 
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Appendix D. Time-varying Ornstein-Uhlenbeck model 


In this section, we discuss how the ideas presented in j|5] and §|6] can be used to computed 
FPT properties for a time varying Ornstein-Uhlenbeck (0-U) model. 

The O-U model captures decision-making through the first passage of trajectories of 
an O-U process (Cox and Miller 1965) through two thresholds. Our calculations for the 


multistage process also help analyze the Ornstein-Uhlenbeck (O-U) model. Similar to the 
multistage process, the n-stage O-U process with piecewise constant parameters is defined 
by 

dx(t) = n(t)dt — A (t)x(t)dt + a(t)dW(t ), x(t 0 ) = x 0l (D.l) 


where 


p(t) = p*, for ti _i <t <U, 
o it) = <7i , for ti -1 <t <U, 

A (t) = A *, for ti -1 <t<U, 

for each i e {1, Due to the extra —A {t)x{t)dt term, the O-U process is a leaky 

integrator, while the original single stage process is a perfect integrator (A(t) = 0). 

Here, leaky integration means that as the noisy signal is integrated in time with ex¬ 
ponentially increasing (A < 0) or decreasing (A > 0) weights on past observations. Such 
exponential weights lead to ‘recency’ or ‘decay’ effects, i.e., the earlier stages (or late stages) 
may have greater influence on the ultimate decision, whereas with the single stage process, 
all of the signal throughout the entire decision period is weighed equally. 


Appendix D.l. The O-U process as a transformation of the Wiener process 

In this section we show how our calculations for the multistage process can be easily 


applied to decision models driven by O-U processes via a transformed Wiener process (Cox 


and Miller, 1965, §5.9). 


Consider the single-stage O-U process 

dx(t) = pdt — A x{t)dt + adW(t), x(t 0 ) = x 0 . 


(D.2) 


The O-U process (|D.2|) can be written as a time-varying location-scale transformation of 

'<7 2 (exp(2At) — 1)' 


the Wiener process (Cox and Miller, 1965, §5.9), i.e., 


x(t) = y( 1 — exp(—At)) + exp(— Xt)xo + exp(— Xt)W 
A 


2A 


(D.3) 


In order to derive (D.3), note that for O-U process (D.2) 


T, 


x(t) = x(0) exp(—At) + ^(1 — exp(—At)) + a / exp(—As)dlU(s), 

^ Jo 

Note that the stochastic process a exp(—As)dW(s) is equivalent to the stochastic process 
crW((l — exp(—2At))/2A) in the sense of distribution. Furthermore, crW((l — exp(—2At))/2A) 
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is equivalent to the stochastic process e _Ai fU(cr 2 (exp(2Af) — 1)/2A) in a similar sense. This 
means that for each realization of the process a exp(—As)dVF(s) there exists an identical 
realization of the process e~ xt W (a 2 (exp(2Xt) — 1)/2A). 

If we define the transformed time by u(t) := ^ _ _i_ l 0 g(l + 

Then 

x 0 + W (u(t)) = (x(t) - 0 exp (At) + j = (x(t) - jj \j 1 + 2X ^ + j. 

We refer to this process as a Wiener process evolving on exponential time scale. 


Appendix D.2. First passage time of the O-U process 

We now consider the first passage time of the O-U process (D.2) with respect to symmetric 
thresholds ±£. If x(t) = ±£, we have x 0 + W{u) = ( ± £ — (() \J 1 + ^ + ((. We denote 

this last quantity by £ ± (u). Consequently, the FPT for x{t) with respect to thresholds ±£ 
is a continuous transformation of the FPT of a Wiener process starting at Xq and evolving 
on the transformed time u with respect to time-varying thresholds at £ ± (u). Since u is a 
monotonically increasing function, the distribution of the first passage time r can be obtained 
from w(r), the FPT distribution of the Wiener process. Furthermore, the lower threshold 
hitting probabilities of the two processes are the same. 

Note that in transforming the FPT problem for the O-U process (D.2) to the FPT problem 
for the Wiener process evolving on exponential time scale, removes all the parameters from 
the underlying process Xq + W{u) and puts them in thresholds £ =t (u) and exponential time 
scale u. In addition to utilizing the results of £j5] to multistage O-U processes, the above 
transformation is also helpful is speeding up Monte-Carlo simulations of the O-U process. 
Since the transformed process evolves on exponential time scale, the Monte-Carlo simulations 
with transformed process should roughly take time that is a logarithmic function of time 
taken by the O-U process (D.2). 

Computation of FPT distributions for the Wiener process with time-varying thresholds 
is, to our knowledge, not analytically tractable. However, the time-varying thresholds can be 
approximated by piecewise constant time-varying thresholds and approximate FPT distribu¬ 
tions can be computed using the multistage model. While the thresholds £ ± are asymmetric 
for the transformed process (i.e. £ + 7 ^ — £ _ ), unlike in the case described for the multi¬ 
stage process; such a case can be easily handled by replacing the expression in (13) and (18) 
with corresponding expressions for asymmetric thresholds (see Douady, 1999; Borodin and 


Salminen 2002). 
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Appendix D.3. Approximate computation of the FPT distribution of the multistage 0-U 
process 


Similar to the transformation described in Appendix D.l 


cess (D.l) for t G [U_i,U) can be written as 


the multistage O-U pro- 


x{t) = ^(1 - exp(—A i(t - ti- 1 ))) + exp(—Aj(t - U_i))x(U_i) 


+ exp(—A,(f - ti -1 


<r?(exp(2Aj(t - ti- 1 )) - 1) 
2A ?; 


(D.4) 


Let Ui(t) = cr 2 (exp(2Aj(f — U-i)) — l)/2Aj. Also, let r t and i, i G {1,..., n} be defined 
similarly to the multistage model. Then, conditioned on a realization of X\_i, the FPT 
problem of the f-th stage O-U process can be equivalently written as the FPT problem of 
the Wiener process W-i + W(ui(t)) with respect to thresholds 





2A iUift') pi 


We can approximate each stage of the O-U process (D.l) by a multistage process representing 


the above Wiener process with time varying thresholds. This sequence of multistage processes 


is itself a larger multistage process that approximates (D.l) and its FPT distribution can 
be computed using the method developed in §[5} Note that this method only yields FPT 
distributions. The expected decision times and probability of hitting a particular threshold 
can be computed by integrating these distributions. 
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